##--------------------------------------##
# a script to calculate mean plot values #

##read in the plots
plots = read.csv("raw exp plot data.csv",as.is=TRUE)
colnames(plots)[2] = "species"

##read in the traits
traits = read.csv("species trait matrix.csv",as.is=TRUE)
colnames(traits)[2] = "species"
colnames(traits)[4] = "growth_form"
colnames(traits)[5] = "lifespan"

###########################################################################
##make a unique list of the plots
##one site at a time - starting with JRG
JRG_plots = subset(plots, plots$Site=="JRG" & plots$treatment%in%c("C","W"))
JRG_plot_names = with(JRG_plots,paste(treatment,block,plotID,Year))
JRG_plot_list = unique(JRG_plot_names)


treatment = NULL
block = NULL
plotID = NULL
Year = NULL
seed_mass = NULL
forb_cover = NULL
grass_cover = NULL
shrub_cover = NULL
annual_cover = NULL
biennial_cover = NULL
perennial_cover = NULL
LL = NULL
LMA = NULL
Nmass = NULL
Narea = NULL
Pmass = NULL
Parea = NULL
Amass = NULL
Aarea = NULL

for (i in 1:length(JRG_plot_list))
	{
	temp = subset(JRG_plots,JRG_plot_names == JRG_plot_list[i])
	temp_with_traits = merge(temp,traits) 

	treatment[i] = temp$treatment[1]
	block[i] = temp$block[1]
	plotID[i] = temp$plotID[1]
	Year[i] = temp$Year[1]
	seed_mass[i] = with(temp_with_traits,weighted.mean(log_seed_mass,Relabund,na.rm=TRUE))
	forb_cover[i] = sum(with(temp_with_traits,Relabund[which(growth_form=="F")]))
	grass_cover[i] = sum(with(temp_with_traits,Relabund[which(growth_form=="G")]))
	shrub_cover[i] = sum(with(temp_with_traits,Relabund[which(growth_form=="S")]))
	annual_cover[i] = sum(with(temp_with_traits,Relabund[which(lifespan=="A")]))
	biennial_cover[i] = sum(with(temp_with_traits,Relabund[which(lifespan=="B")]))
	perennial_cover[i] = sum(with(temp_with_traits,Relabund[which(lifespan=="P")]))
	LL[i] = with(temp_with_traits,weighted.mean(LL,Relabund,na.rm=TRUE))
	LMA[i] = with(temp_with_traits,weighted.mean(LMA,Relabund,na.rm=TRUE))
	Nmass[i] = with(temp_with_traits,weighted.mean(Nmass,Relabund,na.rm=TRUE))
	Narea[i] = with(temp_with_traits,weighted.mean(Narea,Relabund,na.rm=TRUE))
	Pmass[i] = with(temp_with_traits,weighted.mean(Pmass,Relabund,na.rm=TRUE))
	Parea[i] = with(temp_with_traits,weighted.mean(Parea,Relabund,na.rm=TRUE))
	Amass[i] = with(temp_with_traits,weighted.mean(Amass,Relabund,na.rm=TRUE))
	Aarea[i] = with(temp_with_traits,weighted.mean(Aarea,Relabund,na.rm=TRUE))
	}
JRG_plot_means = cbind(treatment,block,plotID,Year,seed_mass,forb_cover,grass_cover,shrub_cover,annual_cover,biennial_cover,perennial_cover,LL,LMA,Nmass,Narea,Pmass,Parea,Amass,Aarea)
site = rep("JRG",nrow(JRG_plot_means))
years_since_start = Year-min(Year)+1
JRG_plot_means = cbind(JRG_plot_means,site,years_since_start)
write.csv(JRG_plot_means, "Jasper Ridge Plot Trait Menas.csv")

##Cool.  That's the matrix of plot means. Now some tests!

##basic tests, with interaction with year
summary(aov(forb_cover~treatment*Year))	##treatment p = 0.5089
summary(aov(annual_cover~treatment*Year))	##treatment p = 0.1586
summary(aov(seed_mass~treatment*Year))	##treatment p = 0.072219
summary(aov(LL~treatment*Year))		##treatment p = 0.330925
summary(aov(LMA~treatment*Year))		##treatment p = 0.07392
summary(aov(Nmass~treatment*Year))		##treatment p = 0.7513988
summary(aov(Narea~treatment*Year))		##treatment p = 0.09145
summary(aov(Amass~treatment*Year))		##treatment p = 0.89543
summary(aov(Aarea~treatment*Year))		##treatment p = 0.989153

##interactions with year and block - more significant!
summary(aov(forb_cover~treatment*Year*block))	##treatment p = 0.400
summary(aov(annual_cover~treatment*Year*block))	##treatment p = 0.0312
summary(aov(seed_mass~treatment*Year*block))	##treatment p = 0.00370
summary(aov(LL~treatment*Year*block))		##treatment p = 0.0787
summary(aov(LMA~treatment*Year*block))		##treatment p = 0.0194
summary(aov(Nmass~treatment*Year*block))		##treatment p = 0.685
summary(aov(Narea~treatment*Year*block))		##treatment p = 0.0167
summary(aov(Amass~treatment*Year*block))		##treatment p = 0.884
summary(aov(Aarea~treatment*Year*block))		##treatment p = 0.986


##################################################################################
##ok, moving on to Konza.  Much of the following code will replicate that from the
##previous section
KNZ_plots = subset(plots, plots$Site=="KNZ" & plots$treatment%in%c("C","W"))

##the Konza data formatting is messed up - for some reason some plots
##labeled as "NoWater" are coded as "W" in treatment
KNZ_plots$treatment = "W"
KNZ_plots$treatment[grep("NoWater",KNZ_plots$plotname)] = "C"

KNZ_plot_names = with(KNZ_plots,paste(treatment,block,Year,Abbreviation))
KNZ_plot_list = unique(KNZ_plot_names)

treatment = NULL
block = NULL
plotID = NULL
Year = NULL
seed_mass = NULL
forb_cover = NULL
grass_cover = NULL
shrub_cover = NULL
annual_cover = NULL
biennial_cover = NULL
perennial_cover = NULL
LL = NULL
LMA = NULL
Nmass = NULL
Narea = NULL
Pmass = NULL
Parea = NULL
Amass = NULL
Aarea = NULL

for (i in 1:length(KNZ_plot_list))
	{
	temp = subset(KNZ_plots,KNZ_plot_names == KNZ_plot_list[i])
	temp_with_traits = merge(temp,traits) 

	treatment[i] = temp$treatment[1]
	block[i] = temp$block[1]
	plotID[i] = temp$plotID[1]
	Year[i] = temp$Year[1]
	seed_mass[i] = with(temp_with_traits,weighted.mean(log_seed_mass,Relabund,na.rm=TRUE))
	forb_cover[i] = sum(with(temp_with_traits,Relabund[which(growth_form=="F")]))
	grass_cover[i] = sum(with(temp_with_traits,Relabund[which(growth_form=="G")]))
	shrub_cover[i] = sum(with(temp_with_traits,Relabund[which(growth_form=="S")]))
	annual_cover[i] = sum(with(temp_with_traits,Relabund[which(lifespan=="A")]))
	biennial_cover[i] = sum(with(temp_with_traits,Relabund[which(lifespan=="B")]))
	perennial_cover[i] = sum(with(temp_with_traits,Relabund[which(lifespan=="P")]))
	LL[i] = with(temp_with_traits,weighted.mean(LL,Relabund,na.rm=TRUE))
	LMA[i] = with(temp_with_traits,weighted.mean(LMA,Relabund,na.rm=TRUE))
	Nmass[i] = with(temp_with_traits,weighted.mean(Nmass,Relabund,na.rm=TRUE))
	Narea[i] = with(temp_with_traits,weighted.mean(Narea,Relabund,na.rm=TRUE))
	Pmass[i] = with(temp_with_traits,weighted.mean(Pmass,Relabund,na.rm=TRUE))
	Parea[i] = with(temp_with_traits,weighted.mean(Parea,Relabund,na.rm=TRUE))
	Amass[i] = with(temp_with_traits,weighted.mean(Amass,Relabund,na.rm=TRUE))
	Aarea[i] = with(temp_with_traits,weighted.mean(Aarea,Relabund,na.rm=TRUE))
	}
KNZ_plot_means = cbind(treatment,block,plotID,Year,seed_mass,forb_cover,grass_cover,shrub_cover,annual_cover,biennial_cover,perennial_cover,LL,LMA,Nmass,Narea,Pmass,Parea,Amass,Aarea)
site = rep("KNZ",nrow(KNZ_plot_means))
years_since_start = Year-min(Year)+1
KNZ_plot_means = cbind(KNZ_plot_means,site,years_since_start)
write.csv(KNZ_plot_means, "Konza Plot Trait Menas.csv")

##basic tests, with interaction with block
summary(aov(shrub_cover~treatment*block))		##treatment p = 
summary(aov(forb_cover~treatment*block))		##treatment p = 
summary(aov(annual_cover~treatment*block))	##treatment p = 
summary(aov(seed_mass~treatment*block))		##treatment p = 
summary(aov(LL~treatment*block))			##treatment p = 
summary(aov(LMA~treatment*block))			##treatment p = 
summary(aov(Nmass~treatment*block))			##treatment p = 
summary(aov(Narea~treatment*block))			##treatment p = 
summary(aov(Amass~treatment*block))			##treatment p = 
summary(aov(Aarea~treatment*block))			##treatment p = 


###############################################################
##Now, Shortgrass Steppe
##For now, treating each plot as a plot, even though the may be blocked
SGS_plots = subset(plots, plots$Site=="SGS" & plots$treatment%in%c("C","W"))
SGS_plot_names = with(SGS_plots,paste(treatment,block,plotID,Year))
SGS_plot_list = unique(SGS_plot_names)


treatment = NULL
block = NULL
plotID = NULL
Year = NULL
seed_mass = NULL
forb_cover = NULL
grass_cover = NULL
shrub_cover = NULL
annual_cover = NULL
biennial_cover = NULL
perennial_cover = NULL
LL = NULL
LMA = NULL
Nmass = NULL
Narea = NULL
Pmass = NULL
Parea = NULL
Amass = NULL
Aarea = NULL

for (i in 1:length(SGS_plot_list))
	{
	temp = subset(SGS_plots,SGS_plot_names == SGS_plot_list[i])
	temp_with_traits = merge(temp,traits)

	treatment[i] = temp$treatment[1]
	block[i] = temp$block[1]
	plotID[i] = temp$plotID[1]
	Year[i] = temp$Year[1]
	seed_mass[i] = with(temp_with_traits,weighted.mean(log_seed_mass,Relabund,na.rm=TRUE))
	forb_cover[i] = sum(with(temp_with_traits,Relabund[which(growth_form=="F")]))
	grass_cover[i] = sum(with(temp_with_traits,Relabund[which(growth_form=="G")]))
	shrub_cover[i] = sum(with(temp_with_traits,Relabund[which(growth_form=="S")]))
	annual_cover[i] = sum(with(temp_with_traits,Relabund[which(lifespan=="A")]))
	biennial_cover[i] = sum(with(temp_with_traits,Relabund[which(lifespan=="B")]))
	perennial_cover[i] = sum(with(temp_with_traits,Relabund[which(lifespan=="P")]))
	LL[i] = with(temp_with_traits,weighted.mean(LL,Relabund,na.rm=TRUE))
	LMA[i] = with(temp_with_traits,weighted.mean(LMA,Relabund,na.rm=TRUE))
	Nmass[i] = with(temp_with_traits,weighted.mean(Nmass,Relabund,na.rm=TRUE))
	Narea[i] = with(temp_with_traits,weighted.mean(Narea,Relabund,na.rm=TRUE))
	Pmass[i] = with(temp_with_traits,weighted.mean(Pmass,Relabund,na.rm=TRUE))
	Parea[i] = with(temp_with_traits,weighted.mean(Parea,Relabund,na.rm=TRUE))
	Amass[i] = with(temp_with_traits,weighted.mean(Amass,Relabund,na.rm=TRUE))
	Aarea[i] = with(temp_with_traits,weighted.mean(Aarea,Relabund,na.rm=TRUE))
	}
SGS_plot_means = cbind(treatment,block,plotID,Year,seed_mass,forb_cover,grass_cover,shrub_cover,annual_cover,biennial_cover,perennial_cover,LL,LMA,Nmass,Narea,Pmass,Parea,Amass,Aarea)
site = rep("SGS",nrow(SGS_plot_means))
years_since_start = Year-min(Year)+1
SGS_plot_means = cbind(SGS_plot_means,site,years_since_start)
write.csv(SGS_plot_means, "Short Grass Steppe Plot Trait Menas.csv")

summary(aov(shrub_cover~treatment*block))		##treatment p =
summary(aov(forb_cover~treatment*block))		##treatment p = 0.001772
summary(aov(annual_cover~treatment*block))	##treatment p = 0.0596
summary(aov(seed_mass~treatment*block))		##treatment p =  0.0567
summary(aov(LL~treatment*block))			##treatment p = 0.00787
summary(aov(LMA~treatment*block))			##treatment p =
summary(aov(Nmass~treatment*block))			##treatment p = 0.08367
summary(aov(Narea~treatment*block))			##treatment p =
summary(aov(Amass~treatment*block))			##treatment p =
summary(aov(Aarea~treatment*block))			##treatment p = 0.01199

plot(as.factor(treatment),Aarea)
plot(as.factor(treatment),LL)


###############################################################
##Now, Sevilleta
##For now, treating each plot as a plot, even though the may be blocked
SEV_plots = subset(plots, plots$Site=="SEV" & plots$treatment%in%c("C","W"))
SEV_plot_names = with(SEV_plots,paste(treatment,block,plotID,Year,Experiment))
SEV_plot_list = unique(SEV_plot_names)

treatment = NULL
block = NULL
plotID = NULL
Year = NULL
seed_mass = NULL
forb_cover = NULL
grass_cover = NULL
shrub_cover = NULL
annual_cover = NULL
biennial_cover = NULL
perennial_cover = NULL
LL = NULL
LMA = NULL
Nmass = NULL
Narea = NULL
Pmass = NULL
Parea = NULL
Amass = NULL
Aarea = NULL

for (i in 1:length(SEV_plot_list))
	{
	temp = subset(SEV_plots,SEV_plot_names == SEV_plot_list[i])
	temp_with_traits = merge(temp,traits)

	treatment[i] = temp$treatment[1]
	block[i] = strsplit(strsplit(SEV_plot_list[i],split = "_")[[1]][1]," ")[[1]][3]
	plotID[i] = temp$plotID[1]
	Year[i] = temp$Year[1]
	seed_mass[i] = with(temp_with_traits,weighted.mean(log_seed_mass,Relabund,na.rm=TRUE))
	forb_cover[i] = sum(with(temp_with_traits,Relabund[which(growth_form=="F")]))
	grass_cover[i] = sum(with(temp_with_traits,Relabund[which(growth_form=="G")]))
	shrub_cover[i] = sum(with(temp_with_traits,Relabund[which(growth_form=="S")]))
	annual_cover[i] = sum(with(temp_with_traits,Relabund[which(lifespan=="A")]))
	biennial_cover[i] = sum(with(temp_with_traits,Relabund[which(lifespan=="B")]))
	perennial_cover[i] = sum(with(temp_with_traits,Relabund[which(lifespan=="P")]))
	LL[i] = with(temp_with_traits,weighted.mean(LL,Relabund,na.rm=TRUE))
	LMA[i] = with(temp_with_traits,weighted.mean(LMA,Relabund,na.rm=TRUE))
	Nmass[i] = with(temp_with_traits,weighted.mean(Nmass,Relabund,na.rm=TRUE))
	Narea[i] = with(temp_with_traits,weighted.mean(Narea,Relabund,na.rm=TRUE))
	Pmass[i] = with(temp_with_traits,weighted.mean(Pmass,Relabund,na.rm=TRUE))
	Parea[i] = with(temp_with_traits,weighted.mean(Parea,Relabund,na.rm=TRUE))
	Amass[i] = with(temp_with_traits,weighted.mean(Amass,Relabund,na.rm=TRUE))
	Aarea[i] = with(temp_with_traits,weighted.mean(Aarea,Relabund,na.rm=TRUE))
	}
SEV_plot_means = cbind(treatment,block,plotID,Year,seed_mass,forb_cover,grass_cover,shrub_cover,annual_cover,biennial_cover,perennial_cover,LL,LMA,Nmass,Narea,Pmass,Parea,Amass,Aarea)
site = rep("SEV",nrow(SEV_plot_means))
years_since_start = Year-min(Year)+1
SEV_plot_means = cbind(SEV_plot_means,site,years_since_start)
write.csv(SEV_plot_means, "Sevilleta Plot Trait Menas.csv")


summary(aov(shrub_cover~treatment*Year))		##treatment p = 0.04115
summary(aov(forb_cover~treatment*Year))		##treatment p = 0.00000000234
summary(aov(annual_cover~treatment*Year))	##treatment p = 0.0000000554
summary(aov(seed_mass~treatment*Year))		##treatment p =
summary(aov(LL~treatment*Year))			##treatment p =
summary(aov(LMA~treatment*Year))			##treatment p = 0.00005584
summary(aov(Nmass~treatment*Year))			##treatment p = 0.001448
summary(aov(Narea~treatment*Year))			##treatment p = 0.00003264
summary(aov(Amass~treatment*Year))			##treatment p = 0.00002755
summary(aov(Aarea~treatment*Year))			##treatment p = 0.0000004696